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ABSTRACT 

We  present  a  new  algorithm  for  solving  two-stage  stochastic 
mixed-integer  programs  (SMIPs)  having  discrete  first-stage 
variables,  and  continuous  or  discrete  second-stage  variables. 
For  a  minimizing  SMIP,  the  BEST  algorithm  (1)  computes 
an  upper  Bound  on  the  optimal  objective  value  (typically 
a  probabilistic  bound),  and  identifies  a  deterministic  lower- 
bounding  function,  (2)  uses  the  bounds  to  Enumerate  a 
set  of  first-stage  solutions  that  contains  an  optimal  solu¬ 
tion  with  pre-specified  confidence,  (3)  for  each  first-stage 
solution.  Simulates  second-stage  operations  by  repeatedly 
sampling  random  parameters  and  solving  the  resulting  model 
instances,  and  (4)  applies  statistical  Tests  (e.g.,  “screening 
procedures”)  to  the  simulated  outcomes  to  identify  a  near- 
optimal  first-stage  solution  with  pre-specified  confidence. 
We  demonstrate  the  algorithm’s  performance  on  a  stochastic 
facility-location  problem. 

1  INTRODUCTION 

Stochastic  mixed-integer  programs  (SMIP)s  arise  in  varied 
contexts  such  as  industrial  capacity  planning  (e.g.,  Stafford 
1997,  Ahmed,  et  al.  2000),  vehicle  routing  or  allocation 
(e.g.,  Frantzeskakis  and  Powell  1990,  Morton  and  Kenyon 
2001),  facility  location  (e.g.,  Laporte  et  al.  1994),  and 
network  interdiction  (Cormican  et  al.  1998,  Israeli  1999). 
The  common  thread  is  that,  in  the  problems’s  first  stage, 
the  user  must  make  discrete,  resource-constrained  decisions 
about  the  configuration  of  a  system,  and  in  the  second 
stage  uncertainty  resolves  itself  and  the  user  (“adversary” 
in  the  case  of  interdiction  problems)  operates  the  configured 
system  optimally.  The  second-stage  decision  variables  may 
be  continuous  or  discrete.  This  paper  standardizes  on  a 
minimizing  SMIP  (e.g.,  minimizing  cost). 

Most  of  the  literature  on  SMIPs  focuses  on  solving 
problems  for  which  all  second-stage  scenarios  can  be  enu¬ 
merated  (Klein  Haneveld  and  van  der  Vlerk  1998,  Ahmed 
2004).  Two  exceptions  include  sequential  approximation 


(SA)  (e.g..  Kail  et  al.  1988)  and  the  sample  average  ap¬ 
proximation  method  (SAAM)  (Mak  et  al.  1999,  Kleywegt 
et  al.  2001).  SA  sequentially  improves  lower  and  upper 
bounds  by  partitioning  the  state  space  of  the  random  pa¬ 
rameters.  Typically,  SA  solves  a  lower-bounding  problem 
defined  across  n  “conditional-average  scenarios,”  and  must 
increase  n  to  refine  the  partition  and  tighten  the  bound. 
But,  the  computational  workload  tends  to  increase  superlin- 
early  in  n.  SAAM  solves  similar,  /(-scenario  problems  but 
with  sampled  scenario  data.  It  too  must  increase  n  for  bet¬ 
ter  accuracy,  and  therefore  suffers  from  the  computational 
difficulities  that  SA  exhibits.  Both  methods  may  require  spe¬ 
cialized  techniques  to  handle  integer  second-stage  variables. 
Our  goal  is  to  develop  an  easy-to-implement  alternative  to 
SA  and  SAAM  that  does  not  suffer  from  their  drawbacks. 

We  develop  a  fundamentally  new  algorithm  for  solving 
SMIPs,  and  call  it  BEST:  Bound,  Enumerate,  Simulate  and 
Test.  BEST  asks  the  user  to  preselect  e  >  0  and  a  >  0,  and 
then  produces  an  e-optimal  solution  with  a  lower  bound  on 
the  confidence  level  of  approximately  1  —  a;  hereafter  we 
refer  to  this  bound  as  the  “approximate  confidence  level.” 
For  most  practical  problems,  the  approximate  confidence 
level  is  likely  to  be  conservative.  Furthermore,  BEST 
may  produce  a  truly  optimal  solution  with  an  approximate 
confidence  level  that  is  strictly  greater  than  1  —  a. 

In  its  basic  form,  BEST  first  computes  a  global  upper 
bound  (typically  probabilistic)  on  the  SMIP’s  optimal  ob¬ 
jective  value  z*.  It  also  identifies  a  lower-bounding  model 
whose  solution,  for  fixed  first-stage  variables,  provides  a 
restricted  lower  bound  on  z*.  BEST  then  applies  these 
bounds  to  enumerate  a  candidate  set  of  first-stage  solu¬ 
tions  that,  with  pre-specified  confidence,  contains  at  least 
one  optimal  solution.  Typically,  the  candidate  set  is  not  a 
singleton,  so  BEST  proceeds  by  simulating  second-stage 
operations  for  each  candidate:  A  Monte  Carlo  simulation 
samples  random  second-stage  parameters,  and  solves  the 
resulting  optimization  problems.  The  algorithm  then  ap¬ 
plies  a  statistical  test — we  use  bootstrapping — to  screen  out 
solutions  that  are  unlikely  to  be  optimal.  If  a  single  can- 
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didate  remains,  that  solution  can  be  declared  to  be  optimal 
with  high  confidence.  Even  if  more  than  a  single  candidate 
remains,  a  second  test  may  enable  us  to  declare  the  “ap¬ 
parent  best  solution,”  i.e.,  the  solution  having  the  smallest 
average  objective  value,  to  be  e-optimal  with  pre-specified 
confidence.  If  that  test  is  not  satisfied,  the  method  specifies 
a  second  round  of  simulation  and  testing  that  makes  such 
a  declaration  valid. 

BEST  places  modest  requirements  on  the  types  of 
SMIPs  it  can  solve.  A  deterministic  lower  bound  must  be 
available  as  a  function  of  the  first-stage  variables.  And,  the 
expected  cost  of  the  SMIP,  given  a  fixed  first-stage  solution, 
should  be  reasonably  easy  to  estimate  using  Monte  Carlo 
simulation.  The  SMIP  should  also  incorporate  “relatively 
complete  recourse,”  defined  below.  Typically,  we  solve 
SMIPs  with  linear  constraints  and  linear  objective  functions, 
but  “linear”  is  not  an  inherent  requirement.  BEST  does  not 
require  the  solution  of  any  multi-scenario  models  as  do  SA 
and  SAAM:  This  is  its  key  computational  advantage. 

We  seek  to  make  the  statistical-testing  portions  of  BEST 
easily  accessible  the  optimization  community,  so  we  pro¬ 
pose  a  novel  vector  bootstrap  approach  to  screen  candidate 
solutions.  This  eliminates  the  need  for  parametric  char¬ 
acterizations  of  the  joint  distribution  of  objective-function 
values  across  candidate  solutions  (e.g.,  constant  variance), 
and  it  means  that  BEST  can  be  implemented  without  spe¬ 
cialized  statistical  functions  or  tables.  In  fact,  it  can  be 
implemented  entirely  within  an  algebraic  modeling  system 
such  as  GAMS  (Brooke  et  al.  1992). 

The  rest  of  the  paper  is  outlined  as  follows.  “Pre¬ 
liminaries,”  Section  2,  specifies  the  general  formulation  of 
SMIP  and  describes  the  stochastic  facility-location  problem 
we  use  throughout  for  illustrative  purposes.  Section  3  out¬ 
lines  the  BEST  algorithm,  and  briefly  describes  the  bounds, 
enumeration  mechanism,  and  statistical  testing  methods  we 
use.  Section  4  presents  computational  results.  Section  5 
provides  conclusions  and  discusses  directions  for  further 
research. 

2  PRELIMINARIES 

2.1  The  Stochastic  Mixed-Integer  Program  (SMIP) 

We  wish  to  solve,  at  least  approximately,  a  two-stage  SMIP 
with  discrete  first-stage  variables 

x  G  X  =  {x  G  Zni  | Ax  =  b,  0  <  x  <  u},  (1) 

and  with  continuous  or  discrete  second-stage  variables 
y  G  Y  =  {y  G  i?”2 1  some  yj  may  be  integer}.  (2) 


Using  tildes  to  identify  random  parameters,  this  SMIP  is 
SMIP  min  Eh(x,  £) ,  where  (3) 

xEX 

_  ~  T 

h(x,  £)  =  c1  x  +  min  f  y  (4) 

ye  Y 

s.t.  Dy  =  d  +  Bx,  (5) 

and  where  £  =  vec(f,  D,  d.  B).  Typically,  u  =  1,  i.e.,  all 
first-stage  variables  are  binary.  We  assume  that,  with  prob¬ 
ability  one,  the  second-stage  problem  in  y  has  a  bounded, 
feasible  solution  for  any  x  G  X.  Thus,  SMIP  is  a  two- 
stage  stochastic  program  with  relatively  complete  recourse 
(Rockafellar  and  Wets  1976). 

We  use  a  simple  example  of  an  SMIP  throughout  the 
paper  to  illustrate  the  BEST  approach,  a  single -product,  ca¬ 
pacitated  stochastic  facility-location  problem  (SFLP).  Birge 
and  Louveaux  (1997,  pp.  57-59)  and  Laporte  et  al.  (1994) 
describe  similar  models: 


Stochastic  Facility-Location  Problem  (SFLP) 

Indices: 

i  £  I  candidate  facilities,  e.g.,  warehouses 
j  G  ./  customer  zones 

Inputs: 

Cj  deterministic  cost  to  construct  facility  i  ($) 
b  maximum  number  of  facilities 
Ui  planned  capacity  of  facility  i  if  built  (tons) 
fij  deterministic  shipping  cost  from  i  to  j  ($/ton) 
dj  random  demand  in  customer  zone  j  (tons) 

Tj  penalty  for  unmet  demand  at  j  ($/ton) 

Decision  Variables: 

Xi  1  if  facility  i  is  built;  0  otherwise  (1st  stage) 
ytJ  tons  shipped  from  i  to  j  (2nd  stage) 

Vj  tons  unmet  demand  at  zone  j  (2nd  stage) 

Formulation 

z*  =  min  E  h(x,  d) 
xE{o,i}in 

s.t. <  b, 
iei 

h(x,  d)  =  y  dXi  + 
iei 

“"lo  51 E  faVii  +  E  riv3 

i£l j£J  j£J 

S.t.  ^  }  yjj  UiX-i, 

jeJ 


where  (6) 

(7) 

VtGl  (8) 
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+  vi  =dji  v  3  €  J  (9) 

iei 

In  the  first  stage  of  SFLP,  we  choose  which  facilities 
to  construct,  but  in  the  face  of  uncertain  future  demands  for 
product.  In  the  second  stage,  actual  demands  are  realized 
and  the  constructed  facilities  ship  to  meet  those  demands 
as  cheaply  as  possible.  A  penalty  is  paid  for  each  unit 
of  unmet  demand.  For  fixed  x,  the  deterministic  version 
of  SFLP  is  a  transportation  problem  with  elastic  demand 
constraints. 

3  THE  “BEST”  ALGORITHM 

BEST  is  based  on  this  self-evident  proposition: 

Proposition  1  Suppose  z"  and  h'  (x)  are  defined 
for  SMIP  such  that  z”  >  z*  with  confidence  1  —  au,  and 
h'(x)  <  Eh(x,£)  for  all  x  £  X.  Enumerate  X  =  {x  £ 
X\h'(x)  <  z”}  and  assume  X  ^  0.  Then,  with  confidence 
1  —  au,  X  contains  at  least  one  optimal  solution  to  SMIP. 
■ 

Our  methods  for  computing  z"  will  never  give  X  ■  0, 
so  we  can  now  provide  a  well-defined  outline  of  BEST. 

Algorithm  BEST 

Input:  Data  for  an  instance  of  SMIP;  confidence  values  au, 
as,  and  a,  for  bounding,  initial  testing  (“screening”),  and 
final  testing,  respectively,  all  chosen  so  that  1  —  (au  +  as  + 
at)  equals  the  desired  overall  confidence  1  —  a;  allowable 
optimality  gap  e  >  0;  initial  sample  size  no- 
Output:  A  solution  x  to  SMIP  that  is  optimal  with  with 
confidence  at  least  (1  —  au)(l  —  cts),  or  is  e-optimal  with 
confidence  at  least  (1  —  au)(l  —  as)(l  —  at). 

{ 

Call  Bound  to  compute  z" ,  an  upper  bound  on  z* 
having  confidence  level  1  —  au\ 

Call  Enumerate  to  find  the  initial  candidate  set  of 
solutions  X  =  {x  £  X\h'(x)  <  z"}.  X  contains  an 
optimal  solution  with  confidence  1  —  au\ 

If  X  =  {x[j] },  set  e  4—  as  <—  at  <—  0  and  go  to  End; 

Call  Simulate  to  generate  samples  n  =  1, . . .  ,no 
of  £,  and  to  evaluate  h(x.  if)  for  each  sample  and  each 
x  G  X; 

Call  Testl  with  observations  from  Simulate  to  screen 
out  convincingly  inferior  solutions,  leaving  the  selected 
subset  X*  C  X.  X*  contains  an  optimal  solution  with 
(approximate)  confidence  at  least  (1  —  au)(  1  —  as); 

If  X*  =  {x^j},  set  e  <—  at  <—  0  and  go  to  End; 

Call  Test2  with  parameter  e  >  0,  input  X*  and  obser¬ 
vations  on  x  £  X*  from  Simulate; 

If  Test2  returns  n+  >  0,  call  Simulate  with  X  "  replac¬ 
ing  X  and  n+  replacing  no,  but  compute  the  apparent 


best  solution  Xm  with  respect  to  all  n+  +  no  observa¬ 
tions; 

End:  Print!  xm,  “is  an”,  e,“-optimal  solution  with 
approximate  confidence”,  (1  —  au)(l  —  cts)(l  —  a*)); 

} 

3.1  Bounds  for  SMIP 

BEST  requires  a  global  upper  bound  on  z*,  and  a  lower- 
bounding  function  on  Eh(x ,  £  )  for  any  x  G  X .  The  number 
of  possibilities  is  large,  and  many  are  problem  dependent, 
so  we  only  discuss  a  few  options  that  apply  to  SFLP. 

3.1.1  Upper  Bounds 

Specialized  deterministic  bounds  could  be  used  here,  for 
example,  the  Edmunson-Madansky  bound  (Edmunson 
1956,  Madansky  1959),  or  a  bound  based  on  dual  restricted 
recourse  (Morton  and  Wood  1999).  We  would  set  au  =  0  if 
such  a  bound  were  used.  However,  the  following,  standard, 
probabilistic  bound  (e.g.,  Mak  et  al.  1999)  applies  to 
most,  if  not  all  SMIPs,  and  is  easily  described  and  computed. 


Procedure  UpperBound 

Input:  Coefficients  and  distribution  parameters  that  define 
SMIP;  sample  size  n„;  confidence  parameter  au. 

Output:  A  probabilistic  upper  bound  on  SMIP  z"  >  z* 
having  confidence  level  1  —  au. 

{ 

Use  a  heuristic  to  identify  a  “good”  first-stage  solution 
x  to  SMIP. 

According  to  the  distribution  of  £,,  generate  nu  random 
samples,  £1;  £2>  ■  •  • ,  and  evaluate  h(x,  for  each; 

Compute  z"  =  iK^kt)  +  tau,nu-iS/^nf 

where  S  is  the  sample  variance  estimator  for  h(x, 
and  tautnu- 1  is  the  upper  1  —  au  quantile  of  the  t 
distribution  with  nu  —  1  degrees  of  freedom; 

Return  z"\ 

} 

The  bound  z"  is  valid  because  i^ft,(x,4)  >  z*  for 
every  x  £  X.  One  heuristic  that  might  be  used  to  obtain  an 
acceptable  x  would  simply  solve  the  expected-value  problem, 
which  is  this  deterministic  MIP: 

min  h(x,  E£).  (10) 

Or,  we  could  solve  an  approximating  problem  with  a  mod¬ 
est  number  of  samples  fe,  /'  =  1, . . . ,  n' ,  taken  from  the 
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distribution  of  £  (e.g.,  Mak  et  al.  1999): 

1  n' 

(id 

e=i 

This  is  also  a  deterministic  MIR  Our  computational  examples 
in  Section  4  exploit  the  latter  technique. 

3.1.2  Lower  Bounds 

Depending  on  where  an  SMIP’s  random  coefficients  appear, 
BEST  may  be  able  to  use  lower  bounds  based  on  Jensen’s 
inequality  (e.g.,  Birge  and  Louveaux  1997,  pg.  140)  or 
dual  restricted  recourse  (Morton  and  Wood  1999).  Because 
h(x,  d)  for  SFLP  is  convex  in  d,  Jensen’s  inequality  applies: 

Proposition  2  For  SFLP,  h'(x)  =  h(x.  Ed)  < 
Eh(x,  d).  I 

Thus,  we  can  compute  a  lower  bound  on  Eh(x ,  d),  for 
any  x,  by  solving  a  single,  deterministic,  elastic  transporta¬ 
tion  problem  with  demands  set  at  expected  values. 

The  need  for  a  good  deterministic  lower  bound  is,  admit¬ 
tedly,  the  most  restrictive  aspect  of  BEST.  Later  in  the  paper, 
we  discuss  how  to  tighten  the  simple,  deterministic  bound 
described  above  for  SFLP.  When  simple  bounds  like  this  do 
not  apply,  we  conjecture  that  probabilistic  lower-bounding 
techniques  will  prove  useful  (Bayraksan  and  Morton  2006). 

3.2  Enumerating  Candidate  Solutions 

Enumerate  in  BEST  requires  that  we  identify  all  first-stage 
solutions  x  to  SMIP  that  satisfy  li'(x)  <  z" .  Given  binary  x, 
the  following  procedure,  which  can  be  implemented  easily 
in  an  algebraic  model  system  like  GAMS  (Brooke  et  al. 
1992),  will  accomplish  the  task. 

Procedure  Enumerate 

Input.  Deterministic  or  probabilistic  global  upper  bound  z" 
for  SMIP;  data  to  define  the  lower-bounding  function  h' (x)\ 
Output'.  Candidate  solution  set  X  =  {x  €  X\ h'(x)  <  z"}\ 
{ 

JL  <—  0; 

Repeatj 

Attempt  to  solve  minxex  h'(x)  for  x; 

If  this  is  infeasible  or  h'(x)  >  z" ,  Return  X; 

X  <—  X  U  {x}; 

Add  a  constraint  to  the  description  of  X  to  eliminate 
x  but  only  x,  e.g., 

y  xi+  y  (i  -  Xi)  <  |/|  - 1 ;  (12) 

i£l\xi=l  iEl\xi—0 

} 


} 

The  problem  minxex  h'  (x)  is  a  sequentially  restricted, 
deterministic,  integer  program  or  mixed-integer  program. 
Clearly,  Enumerate  terminates  finitely. 

The  reader  will  probably  see  that  the  minimizing  prob¬ 
lem  in  Enumerate  need  not  be  solved  exactly.  In  fact,  much 
computational  effort  can  be  saved  by  halting  the  optimiza¬ 
tion  as  soon  as  it  finds  any  x  €  X  satisfying  h'(x)  <  z" . 
Computational  efficiency  may  also  improve  if  equation  (12) 
can  be  replaced  by  a  stronger  constraint.  For  instance,  if  X 
enforces  a  simple  cardinality  requirement,  lTx  =  b.  (12) 
can  be  replaced  by 

xTx  <  6  —  1.  (13) 

If  any  first-stage  variable  is  a  general,  bounded  integer, 
it  can  be  replaced  with  the  standard  expansion  in  terms  of 
binary  variables  (e.g.,  Owen  and  Mehrotra  2002),  and  the 
enumeration  technique  described  above  then  applies.  We 
mention  a  general,  more  computationally  efficient  procedure 
in  section  5. 

3.3  Simulating  Candidate  Solutions 

To  estimate  the  performance  of  solutions  x  £  X  under 
uncertainty,  we  use  common  random  numbers  (CRNs)  to 
simulate  realizations  of  the  random  second-stage  parameters, 
and  then  solve  the  resulting  optimization  models  to  collect 
optimal  objective  values  for  statistical  analysis.  CRNs  result 
in  greater  efficiency  for  statistical  comparisons  when  they 
induce  positive  correlation,  and  we  expect  such  correlation 
for  many  SMIPs.  For  example,  a  pattern  of  generally  high 
demands  in  SFLP  is  likely  to  result  in  high  distribution 
costs  and  unmet  demand  penalties  for  any  set  of  constructed 
facilities. 

Procedure  Simulate 

Input :  Data  to  define  SMIP;  candidate  solution  set  X  having 
confidence  level  1  —  au;  initial  sample  size  no- 
Output :  Randomly  generated  objective-function  data  for 
each  x  £  X. 

{ 

According  to  the  distribution  of  £,,  generate  no  random 
samples,  l  =  1, . . . ,  n0\ 

Evaluate  z/-e  <—  f°r  aU  and  all  i\ 

} 

3.4  Testing  Candidate  Solutions 

There  are  two  steps  to  the  testing  phase.  First,  a  screening 
test,  Testl,  eliminates  solutions  from  X  that  are  unlikely 
to  be  optimal;  we  use  bootstrapping  here.  The  remaining 
candidate  solutions,  the  selected  subset  X*  C  X,  will  con- 
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tain  an  optimal  solution  with  pre-specified  confidence  of  at 
least  (1  —  au)(l  —  aa).  Typically,  \X*\  is  much  smaller 
than  \X\,  and  the  algorithm  can  terminate  immediately  if 
\X*\  =  1.  This  type  of  screening  procedure  is  known  as 
subset  selection  (e.g.,  Bechhofer  et  al.  1995).  The  second 
phase,  Test2,  handles  situations  with  \X*\  >  1,  and  may 
be  optional. 

Procedure  Testl  (Bootstrap  Screen) 

Input :  Candidate  solution  set  X  created  with  confidence 
1  —  au\  objective-function  samples  zm  for  all  k,  i  from 
Simulate;  screening  confidence  parameter  as\  bootstrap 
sample  size  B. 

Output.  Selected  subset  X*  C  X  which  contains  an  optimal 
solution  with  approximate  confidence  level  (1  —  au)  (1  —  ) . 

{ 

Initialize  Wk  0  for  k  =  1, . . . ,  K  =  \X\\ 

For  each  bootstrap  replication  b  =  1, ...  ,B, 

{ 

Generate  no  indices  from  {1, 2, . . . ,  no}  with  re¬ 
placement.  Call  this  set  L\ 

Compute  average,  optimal  objective  values 

z[b)  < - —y^Zki,  k  =  l,...,K;  (14) 

n0 

1(zlL 

Update  the  tally  for  the  “winning”  solution: 

k*  <—  argmin  z^\  Wk *  Wk*  +  1;  (15) 

k  =  l,...,K 


The  selected  subset  X*  has  random  size  that  depends 
on  the  underlying  distribution  of  objective  values.  If  a  few 
solutions  dominate  the  others,  then  £j**|  will  be  small 
even  if  \X\  is  large  and  we  can  use  a  modest  value  of  hq: 
As  in  other  subset  selection  procedures,  the  most  difficult 
situation  from  a  screening  perspective  occurs  when  multiple 
optima  exist.  The  size  of  X*  also  depends  on  B,  n0,  and  the 
confidence  level.  B  >  1000  is  recommended  for  estimating 
percentiles.  We  do  not  yet  know  how  to  choose  no  a  priori, 
but  we  note  that  bootstrapping  in  other  contexts  can  use 
sample  sizes  as  small  as  nine  (Efron  and  Tibshirani  1993). 
Higher  confidence  (i.e.,  lower  as)  increases  the  expected  size 
of  X*.  For  any  particular  problem,  increasing  or  decreasing 
as  over  limited  ranges  may  not  alter  the  identification  of 
solutions  in  X*.  Thus,  the  procedure  is  conservative,  i.e., 
actual  confidence  level  is  typically  higher  than  the  nominal 
one. 

Additional  sampling  in  the  Test2  phase  may  not  be 
required,  even  if  \X*\  >  1.  If  estimated  objective  values 
for  all  x  £  X*  are  sufficiently  close,  a  simple  test  may 
allow  us  to  declare  the  apparent  best  solution  to  be  e- 
optimal  with  required  confidence  1  —  a*  (assuming  that  X* 
contains  an  optimal  solution  to  begin  with).  Test2 — we 
borrow  “Procedure  JCAf”  from  Kim  and  Nelson  (2001) — 
uses  the  initial  simulated  samples  to  determine  whether 
such  a  declaration  can  be  made,  or  if  additional  sampling  is 
required.  In  the  latter  case,  the  procedure  goes  on  to  specify 
a  bound  on  the  number  of  additional  samples  required  so 
that,  when  the  (possibly  new)  apparent  best  solution  is 
identified,  we  can  validly  declare  it  to  be  e-optimal  with 
confidence  1  —  at. 


} 

Let  [1],  [2], ... ,  [I\]  denote  indices  of  candidate  solu¬ 
tions  so  that  turn  >  m)[2]  >  . . ,  >  w\k].  Then 


n0 


fc=i 


1 

n0 


S  —  1 

yi  ^  <  i 

*:= i 


(16) 


Return  X*; 

} 


This  bootstrapping  approach  preserves  the  correlation 
induced  by  CRNs  for  greater  efficiency.  For  x  £  X,  Test2 
provides  a  direct  estimate  of  the  confidence  level  associated 
with  declaring  this  solution  to  be  “best.”  Asymptotically  in 
no,  Wk/B  must  converge  to  0  for  any  non-optimal  solution, 
to  1/m  for  any  of  m  multiple  optimal  solutions,  and  hence 
to  1  for  a  unique  optimal  solution,  if  one  exists.  (The  nature 
of  SFLP  makes  a  unique  optimum  highly  likely.) 


Procedure  Test2  (Select) 

Input.  Selected  subset  X* ,  \X*\  >  1,  having  (approximate) 
confidence  level  (1  —  a„)(l  —  as);  no  objective-function 
samples  from  Simulate  for  each  x  £  X*;  testing  confidence 
parameter  at ;  optimality  tolerance  e  >  0. 

Output.  The  number  of  samples  in  addition  to  no  that  must 
be  applied  to  each  x  £  X*  to  ensure  that  the  apparent 
best  solution  is  e-optimal  with  (approximate)  confidence 
(1  -  au)(l  -  as)(l  -  at). 

{ 

For  all  Xfc,  x^/  £  X* ,  k  ^  k' ,  compute  sample  variances 
for  the  pairwise  difference  using  the  initial  no  samples: 


S2 


kk' 


no 


— j -  y2  {zu  -  zk'e  ~  ( zk ■ 


e=i 


and  set  S2  maxfc^fc/  S/k, ; 

Compute 


Zk'-)) 


(17) 


h2  (no  -  1) 


2at 

1**1  -  1 


-2/(710  —  1) 

- 1 


(18) 
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Compute  the  maximum,  total  sample  size  required: 

.V  —  _/f252/r2J;  (19) 

Return  n+  <—  max{7V  —  no,  0}; 

} 

The  Test2  procedure  is  less  intuitive  than  the  bootstrap 
in  Testl,  but  in  contrast  to  many  alternatives,  it  requires  no 
special  tables  or  calculations.  In  general,  N  will  be  large  if 
|  A’*  |  is  large,  if  the  (positive)  correlation  induced  by  CRNs 
is  low,  or  if  e  is  small.  As  in  Kim  and  Nelson,  we  reuse 
the  initial  data,  but  note  that  the  estimated  savings  in  the 
total  number  of  samples  can  be  substantial  even  without 
this  reuse. 

4  COMPUTATIONAL  RESULTS 

To  demonstrate  BEST’S  empirical  performance,  we  ran¬ 
domly  generate  15  SFLP  test  problems  with  different  char¬ 
acteristics,  solve  them  using  BEST,  and  present  results  in 
Table  1.  These  problems  vary  by  number  of  customer  zones 
and  number  of  potential  facilities.  Customers  and  poten¬ 
tial  facilities  are  randomly  scattered  across  a  rectangle  with 
aspect  ratio  1:3.  Each  shipping  cost  is  proportional  to  the 
Euclidean  distances  a  shipment  must  travel.  Deterministic 
facility-construction  costs,  facility  capacities,  and  penalties 
for  unmet  demand  are  provided  as  inputs,  as  are  expected 
demands  /Zj  in  each  customer  zone  j.  The  actual  demand 
in  zone  j  £  J  is  modeled  as  Unif(/ij  —  fipj .  fij  +  (3pj)  for 
/3  =  0.1,  0.2, 0.4.  We  use  au  =  as  =  at  =  .025  for  all 
problems,  and  desire  a  solution  within  5%  of  the  optimum. 
Upper  bounds  z"  are  computed  in  Bound  using  nu  =  100 
samples  applied  to  a  heuristic  solution  x  computed  by  solv¬ 
ing  (11)  with  n'  =  20.  We  desire  a  relative  optimality  gap 
of  5%,  so  we  set  e  =  0.05z'  =  0.05minxex  h'(x)\ 

As  anticipated,  the  results  obtained  using  CRNs  are 
highly  correlated.  The  average  pairwise  correlations  (not 
shown)  range  from  a  low  of  0.8706  (for  problem  13)  to  a 
high  of  0.9995  (for  problem  1). 


The  easiest  problems  are  clearly  those  with  smaller 
variability  in  the  average  demand.  Indeed,  Bound  and 
Enumerate  yield  only  a  few  candidate  solutions  when  (3  = 
0.1.  But,  even  in  instances  with  many  candidate  solutions, 
Testl  eliminates  all  but  a  handful  of  these.  Eight  of  the  15 
problems  are  solved  completely  after  Testl,  and  only  one 
of  the  remaining  problems  has  more  than  two  candidate 
solutions.  The  estimated  maximum  relative  gap  for  X, 
denoted  Amax,  is  displayed  in  the  table,  and  is  defined  as 
Amax  maxXt ^k,ex\k^k'{zk  -  zk')/z’.  This  may  be 
viewed  as  a  conservative  point  estimate  of  the  true  relative 


gap  associated  with  the  apparent  best  solution,  and  may 
therefore  compared  directly  to  the  desired  maximum  gap 
of  0.05  (i.e.,  5%). 

Column  11  of  Table  1  provides  the  upper  bound,  com¬ 
puted  from  (19),  on  the  additional  sampling  required  if  the 
initial  data  are  reused;  only  problem  15  requires  further 
simulation  and  testing.  (The  value  of  n+  will  tend  to  in¬ 
crease  as  \X*\  increases,  irrespective  of  S,  and  we  note 
that  the  largest  \X*\  occurs  for  problem  15.  However,  n+ 
is  positive  here  primarily  because  S  is  large.)  Column  12 
shows  the  total  time  required  for  the  Bound,  Enumerate, 
and  initial  Simulate  (BES)  steps,  and  Column  13  shows 
the  total  times  required  to  run  the  bootstrap  screening  test 
Testl. 

4.1  Improved  Lower  Bounds 

A  tighter  lower-bounding  function  h'(x)  for  BEST  may  lead 
to  a  reduced  initial  candidate  set  X  and  therefore  reduced 
computational  workload.  The  ideas  of  sequential  approx¬ 
imation  (SA)  can  help  here.  Applied  to  SFLP,  SA  would 
partition  the  state  space  of  the  random  demand  vector  into 
Q  regions  Vq,  q  =  1, ....  Q,  compute  conditional  expecta¬ 
tions  for  each  element,  Aq  =  f?[d|d  £  Vq\,  and  solve  the 
lower-bounding  problem 

Q 

min  h'q  (x)  =  min  ^P{d€P?}ft(x,d,).  (20) 

9=1 

The  value  of  ftg(x)  will  approach  h(x)  from  below  if 
the  partition  is  refined  and  enlarged  appropriately,  i.e.,  as 
Q  increases.  These  optimization  problems  resemble  Q- 
scenario  stochastic  programs,  which,  of  course,  become 
more  difficult  to  solve  as  Q  increases. 

BEST  does  not  need  an  asymptotically  convergent  lower 
bound,  but  a  tighter  one  may  be  useful.  We  find  that 
substantially  tighter  bounds  can  accrue  at  modest  compu¬ 
tational  cost  by  computing  and  exploiting  a  simple  parti¬ 
tion  based  on  quartiles  of  total  demand.  In  particular,  for 
0  =  D0  <  Di  <  D-2  <  D3  <  ZTj  =  oo,  we  define 

Vq  =  {d|£>9_i  <  ^2  dj  <  Dq},  q  =  1,...,4,  (21) 

3<ZJ 

and  then  compute  D\ ,  I)  >  and  D3  so  that  P{d  £  Vq }  = 
0.25  for  q  =  1,...,4.  Defining  Aq ,  q  =  1,...,4,  as 
indicated  above,  and  solving  (20)  gives  the  new  lower- 
bounding  function. 

We  will  not  present  detailed  results,  but  summarize  the 
effect  of  using  the  improved  bound  on  the  most  difficult 
problems,  problems  6,  9,  12  and  15:  Computational  work¬ 
load  for  the  improved  bound  increases  by  at  most  a  factor 
of  two,  but  |  A  |  is  reduced  by  at  least  an  order  of  magnitude, 
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Table  1:  Computational  Results  for  Stochastic  Facility  Location  Problems 


Problem 

Sites 

\I\ 

Zones 

\J  1 

Max 

Facil. 

b 

P 

Bounds 

1*1 

\X*\ 

^max 

n+ 

Computing 

Times^  (secs) 

z" 

BES 

Testl 

1 

10 

20 

5 

.1 

784.7 

801.9 

2 

1 

0 

0 

15.6 

1.3 

2 

.2 

840.0 

5 

1 

0 

0 

17.4 

4.5 

3 

.4 

939.7 

28 

1 

0 

0 

29.7 

30.5 

4 

16 

40 

8 

.1 

873.4 

878.2 

2 

1 

0 

0 

192.5 

1.3 

5 

.2 

916.7 

13 

2 

0.014 

0 

111.3 

13.2 

6 

.4 

1029.0 

173 

2 

0.028 

0 

423.6 

200.0 

7 

16 

40 

9 

.1 

834.3 

856.0 

8 

2 

0.007 

0 

85.3 

7.7 

8 

.2 

885.1 

13 

1 

0 

0 

78.9 

13.3 

9 

.4 

1016.8 

150 

2 

0.015 

0 

324.2 

168.5 

10 

18 

30 

10 

.1 

518.2 

523.1 

2 

1 

0 

0 

22.1 

1.3 

11 

.2 

563.4 

16 

1 

0 

0 

37.3 

16.5 

12 

.4 

687.8 

566 

2 

0.020 

0 

2296.4 

761.7 

13 

20 

50 

10 

.1 

958.8 

966.2 

3 

2 

0.001 

0 

132.8 

1.1 

14 

.2 

1006.3 

72 

1 

0 

0 

255.3 

74.7 

15 

.4 

1155.1 

1201 

5 

0.021 

127 

27458.1 

1673.3 

t  BES,  i.e..  Bound,  Enumerate  and  Simulate,  from  a  1GHz  Pentium  III  computer  operating  under  Windows  2000, 
Testl  from  a  1.5GHz  Apple  PowerPC  G4  operating  under  Mac  OSX. 


so  the  computational  time  for  the  “BES  steps”  of  BEST 
on  each  of  these  problems  reduces  to  at  most  20%  of  the 
original. 

4.2  Reduced  Simulation  Sampling 

The  results  in  Table  1  have  been  produced  using  an  arbitrarily 
chosen  n0  =  50  samples.  But,  bootstrap  applications  often 
involve  only  10-20  samples  (Efron  and  Tibshirani  1993). 
Accordingly,  we  now  explore  BEST’S  behavior  for  no  =  20 
by  rerunning  Testl  using  the  first  20  samples  produced  by 
Simulate  for  each  x  £  X.  Problems  4,  8,  and  15  each  add 
one  solution  to  the  original  set  X* ,  problem  14  yields  two, 
and  no  changes  appear  for  the  others.  Test2  then  requires 
additional  sampling  for  five  problems,  but  the  total  number 
required  is  less  than  45%  of  the  original  in  all  cases.  The 
results  suggest  that  the  total  number  of  simulated  samples 
could  be  reduced  by  over  50%  while  sacrificing  little  in 
accuracy. 

4.3  Simplified  Testing 

A  few  alternatives  exist  do  exist  for  Testl  and/or  Test2.  For 
example.  Nelson  and  Matejcik  (1995)  describe  a  screening 
procedure  that  uses  CRNs,  which  could  be  used  in  lieu  of 
Testl;  see  also  Nelson  et  al.  2001.  However,  we  would 
prefer  to  replace  Test2  with  a  simpler  bootstrap  procedure, 
or  combine  Testl  and  Test2  into  a  single,  simple  bootstrap 
procedure. 


5  CONCLUSIONS  AND  FUTURE  WORK 

We  have  presented  a  new  method  for  solving  two-stage 
stochastic  mixed-integer  programs  (SMIPs);  first-stage  vari¬ 
ables  must  be  discrete,  but  no  conditions  are  placed  on 
second-stage  variables.  The  BEST  algorithm  (Bound,  Enu¬ 
merate,  Simulate  and  Test)  first  uses  bounding  informa¬ 
tion  to  enumerate  a  candidate  set  of  (first-stage)  solutions 
that  contains  an  optimal  solution  with  high  confidence.  It 
then  simulates  the  behavior  of  each  candidate  solution  by 
sampling  random  second-stage  parameters  and  solving  the 
resulting,  simple,  deterministic  problems.  Next,  it  uses 
statistical  tests — we  apply  bootstrapping — to  screen  out  so¬ 
lutions  that  are  unlikely  to  be  optimal.  If  the  screened 
candidate  set  contains  a  single  element,  the  algorithm  ter¬ 
minates.  Otherwise,  additional  sampling  and  testing  may 
be  applied  to  select  a  single  solution  that  is  e-optimal  with 
high  confidence. 

Much  work  remains  to  make  BEST  better.  We  currently 
enumerate  candidate  solutions  by  solving  a  sequence  of 
increasingly  restricted  MIPs.  A  procedure  much  like  branch 
and  bound  would  perform  this  more  efficiently.  The  need  for 
a  deterministic  lower-bounding  function  in  the  enumeration 
step  limits  BEST’S  applicability,  so  a  probabilistic  alternative 
could  prove  useful.  The  bootstrapping  screening  test  proves 
highly  effective  and  is  simple  to  implement.  But,  if  that 
screening  does  not  yield  a  single  candidate  solution,  our 
follow-on  test  is  more  complicated:  Future  research  will 
investigate  simpler  alternatives. 
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